/*------------------------------------------------------------------------
 *  Copyright (C) 2011  Luis M. de la Cruz Salas
 *
 *  This file is part of TUNA
 *
 *  TUNA is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  TUNA is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 ------------------------------------------------------------------------*/

#ifndef _STRUCTUREDMESH_H_
#define _STRUCTUREDMESH_H_

#include <iostream>
#include <string>
#include <fstream>
#include "Meshes/GeneralMesh.hpp"

namespace Tuna {

template <typename T>  class StructuredMesh;
  
/*!
 *  This structure is used to set the precision, that can be \c float,
 *  \c double and \c long double, as well as the dimension of the problem.
 */  
template <typename T>  
struct Typeinfo <StructuredMesh<T> > {
  typedef typename Typeinfo<T>::prec_t prec_t;
  enum { Dim = Typeinfo<T>::Dim };
};

/*! 
 ***************************************************************************
 * General structured meshes on rectangular domains.
 * This is a general class for structured meshes. The template parameter is
 * and adaptor used to control the behavior of the Structured<> objects.
 * For a Finite Volume discretization the extent of the mesh is defined as
 * shown in the next figure:
 * \image html mallas.png "Meshes in 1D, 2D and 3D" width=10cm
 * \image latex mallas.eps "Meshes in 1D, 2D and 3D" width=10cm
 ***************************************************************************
 *  \author  Luis M. de la Cruz [ Sun Oct 28 17:18:28 GMT 2010 ]
 ***************************************************************************
 */
template<typename Tadaptor>
class StructuredMesh : public GeneralMesh<StructuredMesh<Tadaptor> > {

public:
  
/*! Precision for all calculations. */
  typedef typename Typeinfo<Tadaptor>::prec_t prec_t;
/*! Array of ints to define extents. */        
  typedef typename TunaArray<int, Typeinfo<Tadaptor>::Dim >::tiny intTinyArray_t;
/*! Array of floats to define the deltas in each axis. */        
  typedef typename TunaArray<prec_t, Typeinfo<Tadaptor>::Dim >::tiny floatTinyArray_t;
/*! Arrays for numbers.*/
  typedef typename TunaArray<prec_t,1>::huge ScalarField1D;

  StructuredMesh();  
  StructuredMesh(prec_t l1, int e1);
  StructuredMesh(prec_t l1, int e1, prec_t l2, int e2);
  StructuredMesh(prec_t l1, int e1, prec_t l2, int e2, prec_t l3, int e3);
  StructuredMesh(const StructuredMesh& other);
  ~StructuredMesh() {}
    
// The Curiously Recursive Template Pattern (CRTP).
  Tadaptor& asDerived() { return static_cast<Tadaptor&>(*this); }

  inline void resize(prec_t l1, int e1);
  inline void resize(prec_t l1, int e1, prec_t l2, int e2);
  inline void resize(prec_t l1, int e1, prec_t l2, int e2, prec_t l3, int e3);
    
/*!
 *  Returns the dimension of the mesh.
 */
  inline int getDimension() { return this->Dim; }

/*!
 *  Returns the length of the mesh in each axis.
 *  \param i is X for x-axis, Y for y-axis and Z for z-axis.
 */
  inline prec_t getLength(Axis_t i) { return length(i); }

/*!
 *  Returns the deltas of the mesh in each axis.
 *  \param i is X for x-axis, Y for y-axis and Z for z-axis.
 */
  inline prec_t getDelta(Axis_t i) { return delta(i); }

/*!
 *  Returns the nodes of the mesh in each axis.
 *  \param i is X for x-axis, Y for y-axis and Z for z-axis.
 */    
  inline int getExtentNodes(Axis_t i) { return extent_nodes(i); }

/*!
 *  Returns the volumes of the mesh in each axis.
 *  \param i is X for x-axis, Y for y-axis and Z for z-axis.
 */      
  inline int getExtentVolumes(Axis_t i) { return extent_volumes(i); }

/*!
 *  Returns an array with the lengths of the mesh in all axes.
 */    
  inline const floatTinyArray_t& getLength() { return length; }

/*!
 *  Returns an array with the deltas of the mesh in all axes.
 */    
  inline const floatTinyArray_t& getDeltas() { return delta; }  

  inline int getExtentDeltasVols(Axis_t a) { return dxyz[a].size(); }
/*!
 *  Returns the deltas of volumes for the case of non-uniform meshes
 */
  inline ScalarField1D getDeltasVols(Axis_t a) {
    ScalarField1D aux( dxyz[a].size() );
    int i = 0;
    for(viter = dxyz[a].begin(); viter != dxyz[a].end(); ++viter, ++i) {
      aux(i) = *viter;
    }
    return aux;    
  }  

  inline int getExtentDeltasFace(Axis_t a) { return dface[a].size(); }
/*!
 * Returns the deltas of the faces for the case of non-uniform meshes.
 */
  inline ScalarField1D getDeltasFace(Axis_t a) {
    ScalarField1D aux( dface[a].size() );
    int i = 0;
    for(viter = dface[a].begin(); viter != dface[a].end(); ++viter, ++i)
      aux(i) = *viter;
    return aux; 
  }

/*!
 *  Returns an array with the number of nodes of the mesh in each axis.
 */     
  inline const intTinyArray_t& getExtentNodes()   { return extent_nodes; }

/*!
 *  Returns an array with the number of volumes of the mesh in each axis.
 */     
  inline const intTinyArray_t& getExtentVolumes() { return extent_volumes; }

/*!
 *  Function used to calculate the points of the mesh.
 */    
  inline bool calcPoints() { return asDerived().calcPoints(); }

/*!
 *  Function used to calculate the elements of the mesh.
 */ 
  inline bool calcElements() { return asDerived().calcElements(); }  

/*!
 *  Calculate points in edges for non-uniform meshes defined in 
 *  the adaptor classes.
 *  \param a is X for x-axis, Y for y-axis and Z for z-axis.
 *  \param f the function that gives de distribution of the points.
 *  \param p1 the first point of the segment (default: 0).
 *  \param p2 the second point of the segment (default: total length).
 *  \param n number of points on the segment (default: total nodes).
 */
  template<typename Tfunct>
  inline bool calcPoints(Axis_t a, Tfunct f, 
			 prec_t p1 = 0, prec_t p2 = -1, int n = -1)
  {
    if ( p2 == -1 ) p2 = length(a);
    if ( n == -1 ) n = extent_nodes(a);
    return asDerived().calcPoints(a, f, p1, p2, n);
  }  

/*!
 *  Calculate the deltas for non-uniform meshes.
 */
  inline bool calcDeltas(Axis_t a) { return asDerived().calcDeltas(a); }
    
/*!
 *  Make a mesh equal to other existing mesh.
 */
  StructuredMesh<Tadaptor>& operator=(const StructuredMesh<Tadaptor>& other);
    
  bool print();
  bool print(int I, int J = -1, int K = -1);

protected:
  intTinyArray_t extent_nodes;   ///< Extent for nodes in the mesh
  intTinyArray_t extent_volumes; ///< Extent for volumes in the mesh
  floatTinyArray_t length;       ///< Length of the mesh
  floatTinyArray_t delta;        ///< Mesh increments
  //
  // dxyz[X](i-1)   dxyz[X](i)    dxyz[X](i+1)
  //            \         |          /
  //         |<----->|<------->|<--------->|    <--- deltas at volumes   
  //         |       |         |           |
  //             W        P          E                  
  //       --|---*---|----*----|-----*-----|--- <--- the mesh
  //             |   w    |    e     |
  //             |<------>|<-------->|          <--- deltas at faces
  //                 /          \                
  //          dface[X](i)    dface[X](i+1)             
  //                      
  vector<prec_t> dxyz[3];       ///< Array of deltas for nonuniform meshes
  vector<prec_t> dface[3];       ///< Array of deltas for faces in FV method

  typename vector<prec_t>::const_iterator viter;

private:

};

/*!
 *  Default constructor.
 */
template<typename Tadaptor>
StructuredMesh<Tadaptor>::StructuredMesh() {
  for(int i = 0; i < this->Dim; ++i) {
    extent_nodes(i)   = 1;
    extent_volumes(i) = 1;
    delta(i)          = 1.0;
  }
}

/*!
 *  Constructor for 1D.
 *  \param l1 length of the mesh in the x-axis.
 *  \param e1 number of nodes in the x-axis.
 */
template<typename Tadaptor>
StructuredMesh<Tadaptor>::StructuredMesh(prec_t l1, int e1) {
    extent_nodes(0)   = e1;
    extent_volumes(0) = e1 + 1;
    length(0)         = l1;
    delta(0)          = l1 / (e1 - 1);
}

/*!
 *  Constructor for 2D.
 *  \param l1 length of the mesh in the x-axis.
 *  \param e1 number of nodes in the x-axis.
 *  \param l2 length of the mesh in the y-axis.
 *  \param e2 number of nodes in the y-axis.
 */
template<typename Tadaptor>
StructuredMesh<Tadaptor>::StructuredMesh(prec_t l1, int e1, 
					   prec_t l2, int e2) {
    extent_nodes(0)   = e1;
    extent_nodes(1)   = e2;
    extent_volumes(0) = e1 + 1;
    extent_volumes(1) = e2 + 1;
    length(0)         = l1;
    length(1)         = l2;
    delta(0)          = l1 / (e1 - 1);
    delta(1)          = l2 / (e2 - 1);
}

/*!
 *  Constructor for 3D.
 *  \param l1 length of the mesh in the x-axis.
 *  \param e1 number of nodes in the x-axis.
 *  \param l2 length of the mesh in the y-axis.
 *  \param e2 number of nodes in the y-axis.
 *  \param l3 length of the mesh in the z-axis.
 *  \param e3 number of nodes in the z-axis.
 */    
template<typename Tadaptor>
StructuredMesh<Tadaptor>::StructuredMesh(prec_t l1, int e1, 
					 prec_t l2, int e2, 
					 prec_t l3, int e3) {
    extent_nodes(0)   = e1;
    extent_nodes(1)   = e2;
    extent_nodes(2)   = e3;
    extent_volumes(0) = e1 + 1;
    extent_volumes(1) = e2 + 1;
    extent_volumes(2) = e3 + 1;
    length(0)         = l1;
    length(1)         = l2;
    length(2)         = l3;
    delta(0)          = l1 / (e1 - 1);
    delta(1)          = l2 / (e2 - 1);
    delta(2)          = l3 / (e3 - 1);	
}

/*!
 *  Copy constructor.
 */
template<typename Tadaptor>
StructuredMesh<Tadaptor>::StructuredMesh(const StructuredMesh<Tadaptor>& other) 
{
  extent_nodes    = other.extent_nodes;
  extent_volumes  = other.extent_volumes;
  length          = other.length;
  delta           = other.delta;
}

/*!
 *  Resize the mesh in 1D.
 *  \param l1 length of the mesh in the x-axis.
 *  \param e1 number of nodes in the x-axis.
 */
template<typename Tadaptor>
inline void StructuredMesh<Tadaptor>::resize(prec_t l1, int e1) {
    extent_nodes(0)   = e1;
    extent_volumes(0) = e1 + 1;
    length(0)         = l1;
    delta(0)          = l1 / (e1 - 1);
}

/*!
 *  Resize the mesh in 2D.
 *  \param l1 length of the mesh in the y-axis.
 *  \param e1 number of nodes in the y-axis.
 *  \param l2 length of the mesh in the z-axis.
 *  \param e2 number of nodes in the z-axis. 
 */
template<typename Tadaptor>
inline void StructuredMesh<Tadaptor>::resize(prec_t l1, int e1, 
					     prec_t l2, int e2) {
    extent_nodes(0)   = e1;
    extent_nodes(1)   = e2;
    extent_volumes(0) = e1 + 1;
    extent_volumes(1) = e2 + 1;
    length(0)         = l1;
    length(1)         = l2;
    delta(0)          = l1 / (e1 - 1);
    delta(1)          = l2 / (e2 - 1);
}

/*!
 *  Resize the mesh in 3D.
 *  \param l1 length of the mesh in the y-axis.
 *  \param e1 number of nodes in the y-axis.
 *  \param l2 length of the mesh in the z-axis.
 *  \param e2 number of nodes in the z-axis.
 *  \param l3 length of the mesh in the z-axis.
 *  \param e3 number of nodes in the z-axis.
 */
template<typename Tadaptor>
inline void StructuredMesh<Tadaptor>::resize(prec_t l1, int e1, 
					       prec_t l2, int e2, 
					       prec_t l3, int e3) {
    extent_nodes(0)   = e1;
    extent_nodes(1)   = e2;
    extent_nodes(2)   = e3;
    extent_volumes(0) = e1 + 1;
    extent_volumes(1) = e2 + 1;
    extent_volumes(2) = e3 + 1;
    length(0)         = l1;
    length(1)         = l2;
    length(2)         = l3;
    delta(0)          = l1 / (e1 - 1);
    delta(1)          = l2 / (e2 - 1);
    delta(2)          = l3 / (e3 - 1);
}

/*!
 *  Asignation operator.
 */
template<typename Tadaptor>
StructuredMesh<Tadaptor>& StructuredMesh<Tadaptor>::operator=
(const StructuredMesh<Tadaptor>& other) {
  if (this != &other) {
    extent_nodes   = other.extent_nodes;
    extent_volumes = other.extent_volumes;
    delta          = other.delta;
    length         = other.length;
  } else {
    std::cout << "\n +-----+" 
	      << "\n | Tuna::DiscreteGeom::StructuredMesh::operator= : "
	      <<" \n | WARNING : you are trying to do : mesh = mesh"
	      << "\n +-----+";
  }
  return *this;
}

/*!
 *  Print the info of the mesh on to the screen.
 */
template<typename Tadaptor>
bool StructuredMesh<Tadaptor>::print() {

  string axes = "XYZ";

  std::cout << "\n +----- TUNA info -----+"
	    << "\n | StructuredMesh< ";
  asDerived().print();

  for(int i = 0; i < this->Dim; ++i) {
    std::cout << "\n | " << axes[i] << "-axis : length = " << length(i)
	      << "; nodes = "<< extent_nodes(i)
	      << "; volumes = "<< extent_volumes(i) 
	      << "; delta = " << delta(i);
  }
  std::cout << "\n +---------------------+\n ";

  return 0;
}

/*!
 *  Print the info of the mesh on to the screen (parallel version).
 */
template<typename Tadaptor>
bool StructuredMesh<Tadaptor>::print(int I, int J, int K) {

  string axes = "XYZ";

  std::cout << "\n +----- TUNA info -----+"
	    << "\n | StructuredMesh< ";
  asDerived().print();

  if (this->Dim == 1)
    std::cout << "\n | Proc ("<< I << ") ";
  if (this->Dim == 2)  
    std::cout << "\n | Proc ("<< I << " ," << J <<") ";
  if (this->Dim == 3)
    std::cout << "\n | Proc ("<< I << " ," << J << ", " << ") ";


  for(int i = 0; i < this->Dim; ++i) {
    std::cout << " | " << axes[i] << "-axis : length = " << length(i)
	      << "; nodes = "<< extent_nodes(i)
	      << "; volumes = "<< extent_volumes(i) 
	      << "; delta = " << delta(i);
  }
  std::cout << "\n +---------------------+\n ";

  return 0;
}


} // Tuna namespace

#endif // _STRUCTUREDMESH_H
